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Abstract 


The Global Reference Atmospheric Model (GRAM) program includes the 

capability for simulating pseudo-random perturbations in density, temperature, 

pressure, or wind components along a simulated reentry trajectory or other path 

through the atmosphere (e.g. a vertical profile). Some concerns have been 

expressed by GRAM users, however, that mean- square perturbation gradients [e.g. 
2 

<(Af/Az) > for any perturbation parameter f] may be too large for small values 

of vertical separation Az. The present GRAM perturbation simulations, based on 

a one-step autoregressive model [AR(1)], yield a power spectrum versus 

-2 

wavenumber k which is proportional to k at high wavenumbers . This feature 

2 

also produces mean- square perturbation differences [e.g. <(Af) >] which are 
directly proportional to Az, and mean- square perturbation gradients which are 

inversely proportional to Az. Thus, root-mean- square gradients, (Af/'Az) , 

- 1/2 rmS 
increase with decreasing Az as Az . This report suggests simple 

modifications to GRAM (e.g. changes in the form of the autoregressive 

correlation function used) which overcome this problem, i.e. which produce root- 

mean-square gradients that remain bounded as Az approaches zero. Possible 

applications of more sophisticated simulation approaches, e.g. second order 

autoregressive models [AR(2)], or fractal model techniques, were also explored 

briefly but found to yield improvements which appear to be too small to justify 

their considerable added complexity for use in the GRAM program. 
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INTRODUCTION AND BACKGROUND 


Figure 1 illustrates an application of the present GRAM perturbation model 
in simulating vertical density gradients for vertical profiles. Vertical steps 
of 2 km or 0.25 km were used. Larger density gradient (Af/Az) values are 
evident with the 0.25 km spacing. 

The random perturbation model in GRAM is basically a simple Markov process, 
i.e. one which uses a scalar factor, rather than a transition probability 
matrix, to relate the current perturbation value to the previous perturbation 
value. In the notation of time series (Box- Jenkins) models (see, for example, 
Vandaele, 1983), the GRAM perturbation model is a first-order autoregressive 
model, or AR(1) model, with an exponential correlation function p(Ar) = e 
for spatial separation Ar. L is the integral scale of the correlation function 
p(Ar), that is. 


L - p( Ar) d (Ar) . (1) 

Spectra consistent with this correlation function would be of the form 
(Lumley and Panofsky, 1964) 

kF(k)/a 2 - (kL/jr)/(l + k 2 L 2 ) (2) 

for scalar quantities (density, temperature, etc.) and for the longitudinal 
spectra of vector quantities (wind components) . For the transverse spectra of 
vector components, the spectrum would be 

kF(k)/<r 2 - (kL/2w)(l + 3k 2 L 2 )/(l + k 2 L 2 ) 2 . (3) 

-2 

Both of these spectra vary as F(k) a k at large values of k (kL » 1), i.e., 

for small scales of separation. The spectral forms are widely used in 

turbulence simulations, where they are referred to as the Dryden spectrum (see, 

for example, Fichtl, 1977, or Turner and Hill, 1982). 

Figure 2 shows a comparison between the vertical spectrum of horizontal 

wind (equation 3) evaluated from the GRAM perturbation model and spectra 

presented recently by Van Zandt (1985) . For wavelengths less than about 1 km 

- 3 

(wave numbers greater than 10 cycles/m) the observed spectra are consistent 

_3 

with F(k) a k , a characteristic of the spectrum of a saturated field of 
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gravity waves (Smith et al. , 1987). 

The form of the exponential correlation (and its associated Dryden spectral 
form) were originally selected for GRAM from studies of vertical structure 
function (and a limited amount of horizontal structure function) data. Examples 
of these earlier results (Justus and Woodrum, 1972) are shown in Figure 3. For 
vertical displacement, Az , the AR(1) model used in the GRAM perturbation routine 
produces a structure function which is given by (Justus et al., 1986) 

<Af 2 > = < [ f ( z+Az ) - f ( z ) ] 2 > - 2 a 2 [ 1 - p( Az)] , (4) 

where a ^ is the standard deviation of the perturbation values. As shown by 
Figure 3, the structure function model provided by use of the exponential 
correlation [p(Ar) - exp(-Ar/L)] in equation (4) yields a good fit to the 
observed mean-square vertical differences, at least for vertical separations of 
lkm or larger. Structure function values of vertical separations of less than 1 
km are limited by lack of vertical resolution of the Meteorological Rocket 
Network sensors. 

For vertical displacement, Az, the AR(1) model used in the GRAM perturba- 
tion routine produces rms gradients which, from equation (4), are given by 

(Af/Az) rms - a f {2[l - p(Az)] } 1/2 /Az . (5) 

For the exponential correlation function, and Az/L « 1, equation (5) can be 
approximated as 

(Af/Az) rmg - <? f (2/AzL) 1/2 . (6) 

Thus, as the vertical step Az is decreased in the AR(1) simulation, the rms 

perturbation (density or other parameter) gradient will increase inversely as 

the square root of Az . This is the phenomenon illustrated in Figure 1. 

Comparison of the spectra in Figure 2 and the structure function data in 

Figure 3 suggests that, if the vertical spacing is 1 km or greater the 

perturbation results will be consistent with observed spectral magnitudes and 

-2 

gradients. However, the fact that the model spectrum a k versus the observed 
-3 

spectrum oc k for large k (small displacement) suggests that unrealistic 
gradients and spectral magnitudes might result for vertical separation of less 
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than 1 km. For this reason it has been recommended (Justus, et al., 1986) that 
simulations be done using GRAM with a minimum vertical spacing of 1 km. The 
results of the research reported here are intended to improve the perturbation 
simulations with GRAM in such a way as to allow this restriction to be relaxed. 

METHODS EXAMINED FOR REVISED GRAM PERTURBATION SIMULATION 

The correlation function p(Ar) enters into the simulation of any 
perturbation function f(r) [with mean = 0 and standard deviation = via the 
relation for the first-order autoregressive, AR(1) , model (Vandaele, 1983) 

f(r) - p(Ar) f(r-Ar) + a f [l-p 2 (Ar) ] 1//2 a(r) , (7) 

where a(r) is a random variable with mean - 0 and variance - 1, and Ar is the 

step size in successive positions. The AR(1) model produces a statistically 

stationary series of values for f(r) only if p ( n Ar) - p U ( Ar) , that is, if p(Ar) 

is the exponential function p(Ar) - e Ar / L (Vandaele, 1983; Chapter 4). If the 

condition of statistical stationarity is relaxed, by using a correlation 

function other than e then equation (6) will not produce the desired 

standard deviation (a^) for a series of simulated values of f(r) [see, for 

example, Figure 4, which shows the results of using a Gaussian correlation in 

the AR(1) model, and for which the values of the computed standard deviation and 

gradients are too small at the smaller vertical separation of 0.25 km]. 

-2 

However, it is the exponential correlation function which imposes the k large- 

wavenumber spectrum [equations (2) and (3)] and the rms -perturbation gradients 

which become too large as the separation Ar gets small [equation (6)]. 

In order to remedy this situation, one approach examined was the use of a 

modified correlation function which yields bounded values of rms gradients when 

Ar is small, hopefully without significant adverse impact on the statistical 

stationarity and the estimation of the proper variance for the f(r) series. If 

the leading terms in the expansion of p(Ar) are quadratic in Ar, rather than 

linear in Ar as for the exponential correlation, then equation (5) would yield 

rms perturbations which approach a constant value as Ar approaches zero. 

Candidate modifications for the exponential correlation function examined 

2 2 

include the Gaussian [p(Ar) * exp(-7r Ar /4L )], the parabolic [p(Ar) = 1 - 
2 

(2 Ar/3L) for Ar^3L/2; p( Ar) - 0 for Ar>3L/2], and the modified exponential 

2 

[p(Ar) - 1 - a(Ar/L) up to some small value of Ar and p( Ar) = exp(-b Ar/L) for 
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larger values of Ar, with factors a and b determined from continuity of p(Ar) at 
the transition value and the integral scale condition, equation (1)]. Other 
candidate correlation functions examined were sin(7rAx/2L)/(7rAx/2L) and (1 + 
2Ax/L) exp ( - 2 Ax/L) . 

Another approach examined was the suitability of using a two-step 
autoregressive model, AR(2) . For this model, Vandaele (1983) gives the relation 


f (r) = ^fCr-Ar) + ^ 2 f(r-2Ar) + a { [ 1 - 

where p^ is the one -step autocorrelation 
autocorrelation, and the coefficients 


Vi ' V 2 ]1/2a(r) * (8) 

function value, p 2 is the two-step 
and <f > 2 are given by 


^ = p,( 1 - p 0 )/( 1 - Pi) , and 

J- X 

4>i - (P 2 - p[)/( 1 - p\) 


(9) 


The only restrictions on the correlation function p(Ar) [and the one- step and 
two step values, p^ and , it yields] are that, for statistical stationarity , 
it is necessary that 1, 1 , and \* 2 \<i (Vandaele , 1983 ; Chapter 4) . 

By using equations (8), we see that all three of these conditions are met if p„ 

2 2 2 Z 
is greater than p^ or if p ^ satisfies the relation 2p^-l<P2<P^. The AR(2) model 

will thus be stationary for any correlation function which meets these 

conditions. If the correlation function violates the statistical stationarity 

constraint in the AR(2) model, then the term inside the square root in equation 

(8) becomes negative, and no simulations can be calculated by this method. If 

p^ is equal to p^ then the correlation function is exponential [p(Ar)=e Ar / L ] 

and equations (9) reduce to and Thus, the AR(2) model becomes the 

same as the AR(1) model if the autocorrelation function is exponential. 

Another method examined for the perturbation simulations is a fractal 
simulation technique. For the fractal approach, a series of values for 
perturbations f(r) can be simulated (Carpenter, 1980) by 


f (r) - [f (r+Ar) + f(r-Ar)]/2 + (c a Ar/L) a(r) 


( 10 ) 


where c is a constant (the "roughness factor"), and Ar is at first some large 
separation which is then successively reduced in a "recursive midpoint 
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reduction" process. The advantage of the fractal technique is that it can be 
designed to "mimic" natural curves and fluctuations quite closely. The 
disadvantage, in terms of application in GRAM, is that fractal simulation via 
equation (10) requires a "look- ahead" approach (to evaluate the f(r+Ar) values), 
rather than using only one time step in the "past", as in the current AR(1) 
model . 


RESULTS OF THE STUDY 

For study of the various candidate correlation functions, a set of 10 
series simulations of the equivalent of 100 km were done with vertical steps 

between data points of 2, 1, 0.8, 0.6, 0.4 and 0.2 km. For this test a constant 

value of the integral scale L was taken to be 10 km [the GRAM model actually 
uses height -dependent scale values and employs a two-scale model for large-scale 
and small-scale fluctuations (Justus, et al., 1980)]. The AR(1) model of 
equation (6) and the AR(2) model of equation (7) were both evaluated for each of 
the candidate correlation functions. The fractal model of equation (8) was also 
examined for a single selected value of the roughness factor c. 

The correlation models used, introduced above were: 

• the exponential, p{ Ar) = exp(-Ar/L), 

2 

• the modified exponential, p( Ar) - 1 - a(Ar/L) Ar/L <0.05 

p( Ar) - exp(-bAr/L) Ar/L >0.05 

with a - 19.51615854016301 and 
b - 1.00041693941245578, 

2 2 

• the Gaussian, p( Ar) - exp(-?r Ar /4L ), 

2 

• the parabolic, p(Ar) - 1 - (2 Ar/3L) Ar/L <1.5 

p(Ar) - 0 Ar/L >1.5, 

• p( Ar) - sin(kAr)/(kAr) , with k = 7r/2L, and 

• p( Ar) * (l+mAr)exp( -mAr) , with m - 2/L. 

Averages and standard deviations about the averages of the ten trials are 
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shown in Table 1, for the observed standard deviation cr(f) of the simulation 
series, and for the rms gradient value df/dr - (Af/Az) . The input value for 
a(f) was 10% in each case, so the expected average value cr(f) should be 10% for 
simulation series which are statistically stationary. Of course, as discussed 
above, only the exponential correlation p(Ar) = exp(-Ar/L) is expected to yield 
statistically stationary results for the first order AR(1) model. 

The expected values of the rms gradient df/dr, evaluated by equation (5) 
are given for comparison in Table 2. All correlation functions in both the 
AR(1) and AR(2) models are seen to yield df/dr values close to that expected. 
All of the simulation approaches with all candidate correlation functions except 
the modified exponential are seen to yield df/dr values which are considerably 
below that for the present exponential correlation function, even at vertical 
separations larger than 1 km, where the exponential correlation is considered to 
yield accurate results. Only the modified exponential correlation function 
yields values of df/dr which are consistent with the exponential correlation 
results at separations larger than 1 km. The low value of df/dr for the fractal 
simulations could, of course, be increased by using a larger value for the 
roughness factor, c, than the one selected. However, the practical difficulties 
of implementing the fractal approach in the GRAM program, mentioned above, 
essentially rule out its possible application. 

The results of Table 1 also show that: 

• The Gaussian, parabolic, sin(kAr)/(kAr) , and (l+mAr)exp( -mAr) correlation 
functions produce a(f) values which are too low, because of the statistical 
non-stationarity , when used in the AR(1) model. 

• The modified exponential, parabolic, and sin(kAr)/(kAr) correlation functions 
are not suitable for the second order AR(2) model because of difficulties 
with statistical stationarity (dashes in Table 1) . 

• The Gaussian, parabolic and sin(kAr)/(kAr) correlations produce f) values 
which are too low (even when statistical stationarity problems are not 
present) when used in the second order AR(2) model. 

• The cr(f) values are about as good for the (l+mAr)exp( -mAr) correlation in the 
AR(2) model as for the exponential and modified exponential functions. 
However the (l+mAr)exp(-mAr) correlation must be ruled out for its low 
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resultant values of gradient, df/dr. 


• Although the o(f) values for the fractal model results average close to the 
nominal value of 10.0, the computed a(f) values are more variable (larger ± 
standard deviation about the average) than for the exponential correlation. 

• The a (f ) values for the modified exponential correlation are about as good as 
those for the exponential correlation, when used in the AR(1) model, even in 
the separation range Ar/L <0.05 where the correlation function changes to a 
form which no longer preserves strict statistical stationarity . 


• The df/dr values from the modified exponential correlation are consistent 
with those from the exponential correlation for the range Ar/L > 0.05, but 
remain bounded by approaching a constant value for separations smaller than 
Ar — 0.05L, as required. 


With all of these results taken into consideration, only the modified exponen- 
tial correlation function in the AR(1) model meets the necessary requirements of 
consistency with observed a(f) and df/dr values at large separations (Ar > about 
1 km) , while yielding bounded values of df/dr for small separations (Ar < about 
1 km) . The one possible exception would be the fractal model with a larger 
value of the roughness factor, c. However, practical problems of implementation 
in GRAM rule this out without very major program modifications. 

Figure 5 shows simulated vertical profiles of density perturbations 
computed by using the modified exponential correlation function. Comparison of 
Figure 5 with Figure 1, computed with the simple exponential correlation 
function, shows that the gradient values for the 0.25 km vertical separation 
case are significantly reduced, while the rms perturbation values themselves are 
still consistent with the input value. The results of Figures 1 and 5 (as well 
as Figure 4) were all computed by a two-scale perturbation model, similar to 
that actually used in GRAM, with a small scale value of 10 km, and a large scale 
value of 20 km. Input values of rms magnitude for the perturbations were 10% 
for both small scale and large scale perturbations (14.1% for the total 
perturbations shown in the Figures). Figure 5, combined with the results 
summarized in Table 1, thus confirm that the modified exponential correlation 
function adequately addresses the current problems of large gradients at small 
separations, while leaving the perturbation model essentially unchanged (and 
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consistent with observed data) at large separations. 


A STUDY OF HORIZONTAL PERTURBATION STRUCTURE 


During the original study of horizontal and vertical structure functions on 
which the GRAM correlation function and scales are based (Justus and Woodrum, 
1972), only a very limited amount of data was available for horizontal structure 
function analysis. A more extensive data set is that of the density measure- 
ments made on Space Shuttle reentry trajectories. These data were provided to 
us by Mr. John. Findlay of Flight Mechanics and Control, Inc. 

These shuttle data have been used in a structure function analysis by first 
interpolating the data, observed at irregular spacing of horizontal and vertical 
positions, to a constant step size in the horizontal. Since the Shuttle 
trajectories through the 45-95 km range analyzed are nearly horizontal (glide 
slope of about 1/15 or less) , the rms differences between successive values of 
separation are taken to be representative of a horizontal structure function. 

The structure function, given by equation (4) , approaches the value 

<Af 2 > - 2 ct 2 ( Ar/L) , (10) 

for small separations Ar/L, when the exponential correlation function applies. 

The observed Shuttle correlation function data were plotted on log- log scale and 

2 

power- law exponents for <Af > versus Ar were determined by best- fit slopes. 
Data from overlapping height intervals of 45-65km, 55-75 km, 65-85 km and 75-95 
km were considered separately. The average values of the power law exponent 
[expected value 1.0, from equation (6)], were found to be 1.05±0.21, 0.9410.21, 
0.9510.25, and 1.2910.23, respectively, for these four height ranges (1 values 
give the one -standard deviation range about the mean value) . Except for the 
upper height range (where the data are somewhat uncertain) , these results tend 
to confirm the exponential correlation function form for horizontal separations 
as well as for the vertical separations which have been extensively examined 
earlier (Justus and Woodrum, 1972). 
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MODIFICATIONS IN THE GRAM PROGRAM CODE 


The following program code changes will implement the modified correlation 
function model into the GRAM program. Program line number references are as 
given in the GRAM source code listing given in Appendix D of Justus et al. 
(1980). 

After line CORL 45, add the following new function: 


FUNCTION CORREL(X) 

DATA A, B/19 . 51615854016301 , 1 . 00041693941245578/ 
RHO - l./EXP(B*X) 

IF(X .LT. 0.05) RHO - 1. - A*X**2 

CORREL - RHO 

RETURN 

END 


To replace lines PERT 40, PERT 45, PERT 50, PERT 55, PERT 60 AND PERT 65, 
respectively, insert the following new lines of code: 


10 RDS - CORREL (RDS) PERT 40 

30 RTS - CORREL (RTS) PERT 45 

50 RVS - CORREL (RVS) PERT 50 

70 RDL - CORREL (RDL) PERT 55 

90 RTL - CORREL (RTL) PERT 60 

110 RVL - CORREL (RVL) PERT 65 
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MODEL PERTURBATIONS 
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Figure 1 - Simulation of vertical profile of density perturbations 
with the GRAM two-scale model, with two different 
vertical spacings and the exponential correlation function. 
Small scale = 10 km. Large Scale = 20 km. Each with rms 
magnitude, of 10% (total magnitude 14.1%) 
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Figure 2 Comparison of GRAM model wind spectral with observations 
reported by Van Zandt (1985). 
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Figure 3 - Observed vertical structure function of small and large-scale wind and 
temperature perturbations in the 45-65 km height range (Justus and 
Woodrum, 1972). Curves through the data are structure functions computed 
by the exponential correlation function. 
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Figure 4 - Simulation of vertical profile of density 

perturbations with the GRAM two-scale model 
and the Gaussian correlation function, with 
two different vertical spacings. Parameter 
values same as in Figure 1. 
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Figure 5 - Simulation of vertical profile of density 

perturbations with the GRAM two-scale model 
and the modified exponential correlation function, 
with two different vertical spacings. Parameter 
values same as in Figure 1. 
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Table 1 - Computed values of a(f) (nominal = 10.0) and nns value of gradient, df/dr (nominal 
given in Table 2), for first and second order autoregressive models [AR(1) and 
AR(2)], for various correlation functions, and for the fractal model. Dashes 
indicate no results due to violation of the stationarity constraints in the AR(2) 
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Table 2 - Expected values of rms gradient df/dr = (Af/Ar) from 
equation (5) and the various correlation function models 
used in Table 1. 


CORRELATION: 

EXPO- 

NENTIAL 

MODIFIED 

EXP. 

GAUS- 

SIAN 

PARA- 

BOLIC 

sin(kAx)/ 

('kAx') 

(l+mAx)x 
exp( -mAx) 

Ar/L 







0.20 

30.11 

30.11 

12.44 

9.43 

9.05 

17.54 

0.10 

43.63 

43.63 

12.51 

9.43 

9.06 

18.72 

0.08 

49.02 

49.03 

12.52 

9.43 

9.07 

18.97 

0.06 

56.88 

56.89 

12.52 

9.43 

9.07 

19.22 

0.04 

70.01 

62.48 

12.53 

9.43 

9.07 

19.48 

0.02 

99.50 

62.48 

12.53 

9.43 

9.07 

19.74 
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